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A recent reformulation 1] of the problem of Coulomb gases in the presence of a dynam- 
ical dielectric medium showed that finite temperature simulations of such systems can be 
accomplished on the basis of completely local Hamiltonians on a spatial lattice by including 
additional bosonic fields. For large systems, the Monte Carlo algorithm proposed in Ref. [1| 
becomes inefficient due to a low acceptance rate for particle moves in a fixed background 
multiboson field. We show here how this problem can be circumvented by use of a coupled 
particlc-multiboson update procedure that improves acceptance rates on large lattices by or- 
ders of magnitude. The method is tested on a one-component plasma with neutral dielectric 
particles for a variety of system sizes. 

I. INTRODUCTION 

Some years ago, Maggs and collaborators Pi introduced a local reformulation of the statistical 
mechanics of a Coulomb gas of mobile charged entities, relying on the fact that the long-range 
Coulomb interaction could be recovered from a Hamiltonian written in terms of an electric field 
variable coupled to the charge density via Gauss' Law. Since the publication of this paper, more 
recent work has focused on improving the efficiency and capability of the method so that it may be 



applied to larger and more realistic systems 



mm 



|. The partition function in this approach 



is written as a functional integral over an electric field variable containing a physical longitudinal 
(gradient) part incorporating the electrostatic potential energy, as well as an unphysical transverse 
(curl) part. The latter component of the electric field can easily be seen to decouple from the 
charged particle motion, provided the dielectric constant does not change in the course of the simu- 
lation (i.e. is non-dynamical). Note that a spatially varying dielectric medium is perfectly allowable 
in the context of the original method: the integration over the transverse parts of the electric field 
variable only introduces spurious interactions if the dielectric function is also a function of the 
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particle locations. This latter situation is of course quite common in simulations of biophysical 
interest, for example if the charges are attached to a polymer undergoing conformational changes. 

In a recent paper we suggested a new algorithm for compensating for the spurious interac- 
tions introduced by the integration over the curl part of the electric field in situations of this type. 
The extra terms actually correspond simply to the inverse determinant of the Poisson operator 
V • (e(f)V). Removing them therefore amounts to the problem of introducing a positive power of 
the determinant of a local operator into a functional integral. This problem has been intensively 
studied in elementary particle theory, where the inclusion of virtual quark effects in lattice quantum 
chromodynamics simulations require the introduction of just such a determinant in a functional 
integral. In particular, the multiboson method introduced by Liischer jj] seems ideally suited for 
the task at hand, for reasons explained below. Re-expressing the Poisson determinant in terms of 
an integral over a finite number Nb of additional bosonic fields, we arrive at a local Hamiltonian 
and a partition function that can be easily simulated by standard Metropolis/heatbath techniques. 
In this approach, an update of the system proceeds by successive, and independent, updates of 
(a) charged particle locations, (b) electric field values (on the links of the lattice), and (c) all Nb 
multiboson fields (on the sites of the lattice). 

The algorithm studied in Ref. does suffer however from a serious drawback. The multiboson 
fields are in general rather strongly coupled to the particle locations, so that a particle move leaving 
the boson fields unchanged becomes difficult when it induces a large corresponding change in the 
multiboson energy, without allowing the multiboson fields to relax in response to the new particle 
location. The problem becomes acute when the number Nb of multiboson fields is increased: for 
example, to model more accurately the low momentum modes on a large lattice. One then finds an 
acceptance rate for particle moves that decreases exponentially with Nb- The goal of the present 
paper is to show that a modified update procedure which couples particle moves and boson field 
updates greatly ameliorates this acceptance problem. Efficient simulation with a large number of 
multiboson fields is necessary for performing realistic simulations on large lattices. 

In Section 2, we review the basic multiboson algorithm introduced in Ref. valid for systems 
with dynamically varying dielectric functions. The new coupled update method, allowing the 
multiboson fields to relax appropriately as the particles are moved, is explained in detail in Section 
3. Detailed numerical tests of the new method, with comparisons to the original, uncoupled 
approach, are presented in Section 4. Finally, we offer some concluding remarks in Section 5. 
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II. MULTIBOSON FORMULATION OF COULOMB GAS SYSTEMS WITH VARYING 

DIELECTRIC 

In this section we review the multiboson reformulation of Ref. []| for Coulomb gas systems with 
varying dielectric. Consider the partition function of a system consisting of N free charges (mobile 
or fixed) at locations fi, corresponding to a free charge density 

= ^ei5{f- fi) (1) 

i 

where the system is also described by a linear dielectric function (here assumed isotropic) e(r). We 
shall suppress the dependence of e(f) on particle locations, but the procedure we shall describe 
is precisely devised to properly account for this dependence in the electrostatic energetics of the 
problem. Breaking the electric displacement into longitudinal and transverse parts using the general 
Helmholtz decomposition, 

D{f) = -e(f)V0(f) + V X A{f) (2) 
= 511(f) +^*'-(f). (3) 

As the transverse and longitudinal components are orthogonal to each other, one has 

e(f) J """^ e{f) ^ j """^ e{f) ' 

It follows that the imposition of Maxwell's second law (curl-free electric field) implies that the 
electrostatic energy of the system is given purely in terms of the longitudinal part of the electric 
displacement 

Accordingly, the canonical partition function for the system at inverse temperature (5 becomes 

N 

Z= Y[ dfe-^^- (6) 

i=l 

where D^^ must be determined by first solving —V • (eVc/)) = p, from which one obtains -D^ = 
— e(f)V(/)(f). Note that the (non-existent) transverse part of the displacement field plays no role 
in this result. 

In Ref. 3] we have shown how to write a representation for the above partition function by 
writing the Coulomb energy Boltzmann weight e~^^'^^ in terms of an integral over an unrestricted 
electric displacement field D{f) (i.e. one containing both longitudinal and transverse parts), as well 
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as over a set of multiboson fields with a local Hamiltonian. The integration over the multiboson 
fields is arranged to generate a factor which exactly cancels the contribution from the integral 
over the spurious transverse part D^^(r) for arbitrary dielectric functions e(r). Although the factor 
itself is a determinant and highly non-local, the multiboson Hamiltonian is local and therefore 
conventional Monte Carlo simulation techniques can be applied. 

First, we review some notational issues (for further details, see In order to write a well- 
defined functional integral for Z, we consider a lattice Coulomb gas with an electric displacement 
vector field defined on lattice links, where lattice sites are denoted n and /i = 1, 2, 3 indicating 
the spatial direction of the link. Dielectric values are associated with links on the lattice (rather 
than sites) so that the dielectric function on the lattice becomes e^,^. Introducing left (resp. right) 
lattice derivatives A (resp. A), the lattice version of Eq. [21 becomes 

Dn, = D^^^ + D% (7) 

= -en,i^,i4>n + X! ^t^i^'J^y^na (8) 

Here the electrostatic potential is represented by a lattice site field satisfying the Poisson 
equation 

- ^ A^(e„^A^(/)„) = pn (9) 

Pn = ^diSnn (10) 

i 

It is then straightforward to show that the contribution to the functional integral 

Z'= [\[ dDn^e- 2 (11) 

from the integration over the unphysical transverse degrees of freedom yields the factor 
(rin^x/^n/l) det~2(7\/J)j where M. is essentially the lattice Poisson operator: 

6 6 

M\n = (- X! A^Cm^m)-^" = X! f™A„ - ^ eni\n+i (12) 
/i i=l i=l 



The unwanted determinant factor induced by the transverse integrations can be removed by intro- 
ducing a set of additional fields with a local Hamiltonian chosen to generate a positive square-root 
determinant canceling the unwanted contribution. As explained in one begins from a uniform 
polynomial approximation to the function 1/s in the interval [5,V\ for small 6. In terms of the 
complex roots of the polynomial Zk = fJ-k + ii^k, one may choose, for example (the choice is not 
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unique and partly a matter of convenience the Chebyshev polynomial of order 2Nb with 

Nb 

~ Pis) = C 

s 



-I Nb 

- P{s)^Cl[{{s-f^,)' + ul) (13) 

k=l 

1 9'rrk 

= -(l+i)(l-c»^^) (14) 
-k = (15) 

This representation extends as follows to the determinant of a real positive symmetric operator 
with spectrum in the interval [0, 1] 

Nb 

det+5 (7W) ~ n det"^ {{M - fikf + i^l) (16) 

k=l 

where the representation becomes exact in the limit Nb — > oo,6 — > 0. (In practice, the Poisson 
operator needs to be rescaled, so that its spectrum fits in the unit interval.) The crucial step in 
the multiboson procedure is the replacement of the determinant factor in Eq. Elby an equivalent 

(k) 

Gaussian integral over a set of Nb boson site fields (j)n ,k = 1,2, ..Nb- The choice of Nb is dictated 
by the necessity to adequately describe the low eigenvalues of the Poisson operator M: in Ref. Q], 
this was determined empirically by studying a model in which the structure factor of the system 
was analytically known in the charge- free limit. More generally, we expect that as the lowest 
eigenvalue of on a lattice of dimension LxLxL is typically of order 1 /L^ , the value of Nb will 
have to increase as the lattice size does. However, the precise character of the scaling of Nb with 
L clearly depends on the sensitivity of the observables studied to the infrared spectrum of M. 

Inserting a delta- function to enforce the lattice version of Gauss' Law (thereby fixing the lon- 
gitudinal part of the electric displacement appropriately), the expression found for the correct 

n 

partition function (Eq. (1|] becomes 



Z = UdnU dDnf, n #n ^ '5(E ^M^nM " Pn) 

Note that the effective Hamiltonian appearing in the exponent here is completely local. Monte 
Carlo simulation of the system is in principle straightforward: a state of the system is characterized 
by (i) particle locations fi, which in general will influence the dielectric function e(r) (which more 
accurately should be denoted e(r;rj)), (ii) the lattice displacement field Dn^ (respecting locally 
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the Gauss' Law constraint), and (iii) the Nb auxihary scalar fields and simulation requires 
updates of all of these variables with a procedure respecting detailed balance and according to the 
indicated Boltzmann weight. 



III. EFFICIENT MONTE CARLO PROCEDURE FOR COULOMB GAS/MULTIBOSON 

SYSTEM 



In m, the multiboson approach was tested on a system of particles (both neutral and charged) 



with a dielectric different from that of the ambient medium. This model, first studied in Ref. 
provides a clean testbed for algorithms dealing with the full electrostatic energetics of a system 
with a dynamical dielectric function. The latter changes through the simulation as the dielectric 
constant along a link in the direction from a lattice site n is defined through the relation 

^ = 1 + ^. (18) 

where and e„+^ are either the background dielectric constant or the particle dielectric constant 
depending on whether there is a particle on site n or site n + ^ respectively. For example, in the 
simulations discussed below, the background dielectric constant is set to 1.0, while the dielectric 
constant associated with the particle is set to 0.05. If the particles are electrically neutral, then 
we expect the electric displacement to vanish irrespective of the particle locations, leading to a flat 
(wave-number independent) density-density structure factor. The simulations of Eq. Elin Ref. P| 
show that with a relatively small number of multiboson fields {Nb=^ on a 16x16x16 lattice) the 
neutral system yields a flat structure factor within statistical errors. For the rest of this paper, we 
shall use this model system to discuss issues of acceptance rates and simulation efficiency in the 
multiboson approach. 

Unfortunately, the simulation procedure of Ref. P|, which we now review briefly, leads to an 
acceptance rate for particle moves that falls extremely rapidly as the number of multiboson fields 
is increased. For example, in our 16x16x16 test lattice the procedure of Ref. Q gives a particle 
move acceptance rate of 0.011 when 8 multiboson fields are used, but the acceptance rate drops to 
less than 2e-7 when 32 multiboson fields are used. As we emphasized above, increasing Nb will be 
necessary if we wish to go to much larger lattices, in order to apply the method to systems of real 
biophysical interest, for example. Specifically, the Monte Carlo procedure of consisted of 

1. A Metropolis move of a particle (randomly chosen) to a neighboring site, holding the electric 
displacement field Dnfi and the multiboson fields (jxh"* fixed. 
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2. A heat-bath (or worm j^) update of the electric displacement field Dn/j., with particle loca- 
tions and multiboson fields held fixed. 

3. A Metropolis (or heat-bath) update of the multiboson fields, with electric displacement field 
and particle locations held fixed. 

The acceptance problems of the method arise from the first step above: Metropolis moves of 
the particles along lattice links. As the number of multiboson fields increases, the change in the 
contribution of the multiboson Hamiltonian to the Boltzmann weight increases (roughly linearly 
in Nb) and the acceptance rate for moves decreases exponentially. The origin of the decrease 
is simply that the displacement of a particle to a neighboring site leaving the multiboson fields 
unchanged places the system on the tail of the old multiboson Gaussian, with the multiboson 
Gaussian corresponding to the new particle location having very little overlap with the old one. The 
heat-bath approaches used in the second and third steps above have of course no such acceptance 
problems. 

The displacement of the multiboson Gaussian induced by a particle move holds the key to 
finding a more efficient move algorithm: one clearly needs to update the particle positions and 
multiboson fields simultaneously. We now show that the Gaussian dependence of the Hamiltonian 
on the multiboson fields allows us to do this with a heat-bath approach that yields greatly increased 
(for larger values of Nb, by several orders of magnitude) acceptance rates for particle moves. Let us 
denote the multiboson contribution to the Boltzmann weight in Eq. (|17j) by e~^'^^°^. Now consider a 
particle move across a link connecting two nearest neighbor sites a, b on the lattice (i.e. we assume 
that initially there is a particle at one, but not both of these sites). The complete set of multiboson 
fields on the lattice can be grouped as ((/)^, cpl, cp' *), i = 1,2, ..Nb, where the fields not on sites a or 
b are collectively labeled (p'. The next step is to expose the dependence of the multiboson energy on 
just the fields at sites a and b: this will allow us to derive a heat-bath procedure for simultaneously 
moving the particle between a and b and updating all multiboson fields on these two sites, thereby 
"adapting" the multiboson fields dynamically to the updated particle position. As the multiboson 
energy is a quadratic function of the multiboson fields, we must have 
Nb 

Smho. = Y^iAi^i? + Bi{4f + 2Gi</>>| + Ci0^ + Ac/'l) + Snon-ab(</'', e) (19) 
i=l 

Here, the coefficients Ai, Bi,Gi,Ci, Di, as well as the residual term 5non-ab all depend on the 
fields (p' not on the pair of sites a, b, as well as the dielectric field e, which of course depends on 
particle locations via Eq. ()18p . It is straightforward to find explicit expressions for these coefficients 
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in terms of e and the ^' fields. We shall use Greek indices p,a,T,.., running over six values, to 
denote positive and negative spatial directions out of a given lattice site, while unit vectors in the 
corresponding directions are denoted p etc. Now let us assume that site b is the nearest neighbor 
of site a in the a direction. Defining the site field 

En = J2^np (20) 

p 

one finds for the quadratic coefficients 

A, = {E,-^,^)^ + Y^el^ + uj (21) 
p 

Bi = {E,-fiif + Y^el + uf (22) 
p 

Gi = -eaa{Ea + Ei-2lli) (23) 
while for the coefficient linear in one has 

Ci = -2^{Ea + Ea+p-2lli)eap(t>\j^p + 2 ^ ^apea+p,T4>\+p+f (24) 
p^a p+T7^0 

with a similar equation (replacing a by h) for the coefficient Dj. 

Defining diagonal matrices Aij = AiSij, Bij = BiSij,Oij = GiSij, and a {2Nb):>^{2Nb) matrix 
M by 




and a 2Nb dimensional column vector containing the coefficients Ci,Di 




with a similar representation for the fields ^ai4'b' 



$ = 



06 

the Boltzmann weight for the fields 0a > 06 may be written 

g-*M<i>-v<i> ^ g-(*+|yM-i)M($+iM-ii/)+i\/M-iy ^25) 
= e-™ei^^"'^ (26) 
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with 

^ = $ + ^Af-V (27) 

As the multiboson fields for different index i do not interact, the {2Nb)x(2Nb) matrix algebra 
implicit in the above equations actually decouples into Nb independent 2x2 problems. In particular, 
the total weight for the (pa, 4>b field pair, with a particle at either a or 6, is found by integrating 
out these fields: 

Wab = J d$e-**^*-^* (28) 
= det-5(M)e3^^"'^ (29) 

i 

With these preliminaries, we can now state the procedure needed to implement a coupled heat- 
bath particle move/multiboson field update. The first step is to determine a relative Boltzmann 
weight for placing a particle either at site a or site b. This weight is determined as a product of 
three factors: 

1. If the particle is charged, a contribution from the electric field term e ^ ^"'^ where 
the prime indicates that only lattice links connected to site a or site b are included. Here 
the dependence on particle location is entirely through the change a particle move induces 
in the dielectric field. 

2. A contribution from the dielectric factor e 2 Z^n^ os(<^nfj^)^ with the same interpretation of 
the primed sum. 

3. A contribution from the multiboson energy VFa^e"'^"™-'''' arising from the total integrated 
weight Wab of the (pa, (pb fields computed in Ea. ipH)) . together with the portion of the multi- 
boson energy depending on the neighboring multiboson fields (p' (which also depends on 
particle location only through the dielectric field). 

A decision on a final location for the particle (either site a or site b) is now made on the basis 
of the relative weight determined from the product of the three factors described above. Once a 
particle location has been decided, a new set of multiboson fields on sites a and b are calculated 
by heat-bath by determining ^' according to the Gaussian Boltzmann weight e~*^^*, and then 
setting 

/ . 
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FIG. 1: Structure factor for neutral system calculated using the coupled algorithm. Results are shown for 
4, 8, and 16 bosons. 160,000, 400,000 and 800,000 measurement sweeps were used for the runs with 4, 8, 
and 16 bosons respectively. For clarity, only every fifth wavenumber q is plotted. 



IV. SIMULATION RESULTS FOR THE COUPLED HEAT-BATH METHOD 



We tested the efficiency of our improved simulation technique by simulating a system of 1000 
neutral particles on a 16x16x16 lattice as studied in Ref. 0. For the simulations of this paper 
we have chosen a dielectric constant ratio between the particles and the background medium to 
be 0.05 (in typical molecular dynamics biophysical simulations, for example, the protein interior 
is given an effective dielectric constant of ~ 4 while the aqueous medium is ~ 80). In particular, 
in the simulations the particles have a dielectric constant of 0.05 and the background dielectric 
constant is 1.0. The dimensionless inverse temperature is 0.25, as in the previous simulations of 
this system. 

The relevant physical quantity for systems of this type is the structure function defined as the 
Fourier transform S{q) of the coordinate space density-density correlation function (Ref. ^). As 
we work on a spatial lattice, the allowed momentum magnitude values g = |g| are discrete: we have 
chosen to average the structure factor S{q) over degenerate values of the magnitude of the Fourier 
momentum vector (on a 16x16x16 lattice, there are 115 nonzero distinct allowed values). This data 
set is measured every 10 sweeps, and error bars extracted at the end of the runs by blocking the 
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FIG. 2: Acceptance rates for particle moves as a function of number of multibosons. Shown for both earlier 
algorithm where the multiboson updates and particle updates are uncoupled and the improved coupled 
algorithm. 

measurements in order to take autocorrelations into effect. This allows a measurement of a chi- 
squared goodness-of-fit of the final averaged data to the exact answer, which for neutral dielectric 
particles is simply a flat structure function 

where is the number of particles and V the number of spatial points on the lattice. 

As emphasized in Section 2, the choice of the number Nb of multiboson fields introduced to 
model the dielectric contribution to the partition function is dictated by the need to properly 
describe the low (infrared) part of the spectrum of the Poisson operator M, which contains eigen- 
values of order 1/L^ on LxLxL lattices. Too small a choice for Nb can be expected to lead to 
systematic deviations in the structure factor of Eq.|^for small momenta q. On the other hand, too 
large a choice leads to rapidly falling acceptance rates for multiboson updates and is clearly more 
computationally costly per Monte Carlo update step. The systematic deviation in the structure 
factor resulting from too small a choice for Nb is shown in Fig. ^ where we compare the structure 
factors on a 16x16x16 lattice for Nb =4,8,16. It is clear that for this problem, 4 multiboson fields 
are inadequate, giving a distinct curvature to the structure factor, and a large chi-squared per 
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degree of freedom of 36.4, whereas the chi-squared per d.o.f. for the simulation with 16 multiboson 
fields is only 1.94. 

Depending on the physical parameter regime, and the size of lattice used, it will evidently 
be necessary to increase Nb to achieve adequate accuracy, at least for long-range features of the 
physics. Unfortunately, the update procedure introduced in Ref. suffers from the drawback 
that, for fixed physical parameters (lattice size, dielectric constant, etc), the acceptance rates for 
multiboson updates decreases exponentially rapidly with Nb- This makes the uncoupled multi- 
boson method impractical for much larger lattices than those studied in our previous work. In 
Fig. 121 we compare the acceptance rates using the original update procedure, in which particle 
moves were decoupled from updates of the multiboson fields, with the acceptance rates using the 
new procedure described in Section 3, in which particle moves are coupled to readjustments in the 
multiboson field. For larger values of Nb it is clear that the acceptance rates are increased by 
orders of magnitude (e.g. by about 10^ for Nb=32). 

Although the new, coupled procedure is more computationally intensive, it does not require 
significantly more computational time per Monte Carlo sweep. We observe that for a run with 
the parameters described above with 16 bosons, the computational time per Monte Carlo sweep 
using the coupled update method is only about 13% greater than in the original uncoupled algo- 
rithm. This additional computational cost is insignificant given the great increase in particle move 
acceptances that the coupled method provides. 



V. CONCLUSION 



By using an update method where the proposed particle moves include a relaxation of the 
multiboson field to adjust to the change in particle location we have been able to dramatically 
increase the acceptance rates for larger systems that require a larger number Nb of multiboson fields 
to accurately model the long wavelength physics of the system. For both the original uncoupled 
and the improved coupled update algorithms, and for the physical parameters used here, the 
dependence of the acceptance rate on Nb is exponential {exp {—KNb))- The coefficient in the 
exponent, K, is ~ 0.13 for the coupled algorithm and ~ 0.46 for the original method, amounting to 
a difference in acceptance rates of 4 orders of magnitude for a simulation with 32 multiboson fields. 
Tests done with charged particles show similar improvements to the acceptance rates, although for 
some parameter regions the electric field must also be coupled into the particle moves to obtain 
reasonable acceptance rates 
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